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Abstract 

JO The coupled map lattice by Olami et al. [Phys. Rev. Lett. 68, 1244 (1992)] 

^| is "doped" by letting just one site have a threshold, T^ ax , bigger than the 

c5 others. On an L X L lattice with periodic boundary conditions this leads to 

a transition from avalanche sizes of about one to exactly L 2 , and after each 

ON 

avalanche stresses distributes among only five distinct values, r k , related to 
the parameters a and T^ ax by T k = kaT^ ax where k = 0,1,2,3,4. This result 
is independent of lattice size. The transient times are inversely proportional 

o 

to the amount of doping and increase linearly with L. 
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During the last few years there has been a great amount of interest in a coupled map 
lattice introduced by Olami et al. [1]. The map can be derived from a spring-block model 
of an earthquake fault proposed by Burridge and Knopoff [2]. The model consists of blocks 
coupled to its nearest neighbours by elastic springs. Each block is exposed to a uniformly 
increasing stress which discharges when it reaches a certain threshold, and part of the stress 
is then transferred to the nearest neighbours through the couplings. The corresponding 
coupled map lattice is defined on an L X L lattice. Each site is associated with a 

stress T{j. Initially, the stresses are randomly distributed over the lattice. The stresses are 
increased at a slow rate, the same rate for all sites. Finally, one site will reach the threshold, 
T*, and topple. The following relaxation rules are then applied to the system: 

-^nn ^ -^nn ~\~ QL^ij (l) 

T - — >■ 

where nn denotes nearest neighbours to the toppling site The parameter a controls 

the level of conservation in the system. One must have < a < 0.25. The redistribution of 
stress may cause some of the nn sites to exceed the threshold and thus start an avalanche 
of toppling sites. 

For open boundary conditions the model displays a power law in the distribution of 
avalanche sizes for 0.05 < a < 0.25 [3]. The system seems to organize itself into a critical 
state without having to fine-tune the parameter a, and it is therefore claimed to exhibit self- 
organized criticality [4]. This critical behaviour is dependent on boundary conditions and 
on the degeneracy of the system [5,6]. With periodic boundary conditions, the degenerate 
system goes into a state with average avalanche size only slightly bigger than one for almost 
all initial conditions, and it seems clear that the spatial inhomogeneity caused by open 
boundary conditions is responsible for the scaling in this model. The effect of the boundary 
has recently been studied in detail [7]. It has been conjectured that similar avalanche scaling 
should be observed also with other sources of spatial inhomogeneity [5] . 

In order to remove the degeneracy (and thereby introduce a spatial inhomogeneity) we 
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follow Janosi and Kertesz [8] and let the thresholds be site dependent. Site dependent 
thresholds have also been used by Rundle and Klein [10] in a model somewhat similar to the 
model by Olami et al. However, instead of a uniform distribution we use a gaussian distribu- 
tion with mean value 1.0 and standard deviation a. This extension of the model completely 
changes its behaviour both for periodic and open boundaries. For periodic boundaries and 
for several values of a ~ 10~ 2 , our numerical simulations show that the system enters into a 
state with period one and avalanche size exactly equal to L 2 , that is, exactly all sites topple 
once in each avalanche. This is a behavior similar to that in models of globally coupled 
biological oscillators [9]. Furthermore, for a small, the distribution of stresses right after an 
avalanche exhibits only five distinct values, 7% = kaT^ ax where k = 0, 1,2,3,4 and T^ ax is 
the highest threshold in the system (Fig. 1). Only one site holds the largest value, r 5 , and 
this is precisely the site which has the highest threshold. Consequently, this is the site which 
triggers the next avalanche. Thus, all transport of stress in the system must be in the form 
of "packets" of size aT^ ax , and all sites have the same stress T^ ax when they topple. We have 
not observed any case when this did not happen for lattice sizes in the range 15 < L < 50, 
0.05 < a < 0.249 and for many different initial conditions. The tendency towards clustering 
around four values has already been observed by Grassberger, in the limit a — > 0.25 of the 
degenerate model [5] and by Zhang in a slightly different model [11]. Recently, a similar 
"quantization" of stress was observed by Corral et al. in a nonlinearly driven model with 
open boundary conditions [12]. 

The site with the largest threshold is able to trigger all the avalanches only if the difference 
between T^ ax and the stress at this site is less than the difference between threshold stress 
and stress at all other sites. The next highest stress value occuring in the ordered state is 
r 3 = 3aT^ ax . Thus, the relation 
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where T^ in is the smallest threshold in the lattice, seems to be a sufficient condition for the 
strongly ordered period one state, and this is consistent with our numerical results. The 
necessary condition will depend on initial conditions. When a is increased such that this 
condition is violated, the system enters a region of more complicated states, including peri- 
odic behavior with periods larger than one and states with either a very large or an infinite 
period. In the latter case, the avalanche distribution function is exponentially decreasing. 
There seems to be no upper limit on a, other than the conservative limit, a < 0.25, and no 
lower limit on a for obtaining period one. 

One may ask, how many period one attractors of the type described above exist for a 
given lattice? We cannot answer that question, but it is certainly a large number. We 
simulated a lattice with L = 4 and a = 0.01 and 200 different initial conditions ended in 
200 different final configurations — all of the type described above. 

With open boundary conditions, Janosi and Kertesz found an exponential decay of the 
distribution function for a large, uniform spread on the thresholds [8]. For small <r, we 
still find a strong tendency towards very large avalanches. In this case, the distribution of 
stresses immediately after a big avalanche is very similar to the case with periodic boundary 
conditions, but instead of just five values there are now extremely narrow distributions 
around the four lowest 7% values, and in addition there are a few values scattered in the 
intervals between the 7% values (Fig. 2). On inspection, it turns out that the sites with these 
values are all on or near the boundary. 

In order to demonstrate the effect more clearly, we use periodic boundary conditions 
and let just one site have a threshold, T^ ax , larger than the threshold T* of the remaining 
sites. The only difference from the behaviour described above was somewhat longer transient 
times. The time dependence of the avalanche size is illustrated in Fig. 3. The period one 
state implies that the doped site has the stress value r 4 = 4aT^ ax right after an avalanche, 
and that the stresses on the other sites are distributed among the four remaining T^-values. 
It also implies that the doped site triggers all the avalanches. We let T* = 1, and must then 
have 
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(4) 

or 

1 < ^ ax < (5) 
1 — a 

in order to obtain the period one state. For instance, for a = 0.1825 we must have 
1 < T,* < 1.223... . For this value of a we obtained a period one state for T,* < 1.21, 

max j. max — > 

but not for any T^ ax > 1.22 (Fig. 4). In the latter case we observe, as for large <r, a variety 
of different states for different parameter values and initial conditions. In some cases, the 
avalanche distribution function can be fitted to a finite-size scaling hypothesis, but we have 
not been able to find this behavior in any range of parameter values. 

When inequality (5) is fulfilled, it is implied that the distribution of thresholds in the 
interval ((1— ct)T*„„, Tl„) is irrelevant as far as the character of the distribution of thresholds 

\ \ / max ' max/ 

are concerned. However, it may have some influence on the number of sites having each of 
the four possible values. 

We have mainly used lattice sizes up to 50x50, which was also the largest used by 
Olami et al. [1] in their original paper. It has been pointed out [5,13] that the conclusions 
of Olami et al. concerning the finite-size scaling of the avalanche size distribution cannot 
remain correct for sufficiently large lattices. However, in the present case it is observed 
that for periodic boundary conditions, the essential effect of a strongly ordered final state is 
completely independent of lattice size (except, of course, that transient times increase with 
lattice size). Thus, there is no reason to suspect that larger lattices will behave differently. 
However, we have made a single run with a lattice size 200 X 200 and indeed found the same 
strong ordering effect. 

The transient times increases as T,* — > 1 from above. For a = 0.05, 0.10, 0.15 and 0.20, 

max I/ i 

and for T^ ax just above 1, we find that the transient times obeys a power law 

ftr ~ (C ax - I)"", (6) 

where fi = 1.00 ±0.01 (Fig. 5). For a lattice size 15x15, fi is independent of a within errors, 
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and for a = 0.20 it is also independent of L for L = 15,25 and 35 (Fig. 6). The transient 
times seems to increase linearly with L (Fig. 7); the dependence on a is not that simple 
(Fig. 8). 

A strongly ordered state at one value of T^ ax can be used to find a strongly ordered 
state at another value of T^ ax by simply shifting the stresses kaT^ ax in accordance with 
the change in T^ax- According to Eq.( 6) the transition time to a strongly ordered state for 
T'max = 1 (i.e. the original model) would be infinite, but using initial conditions constructed 
in the way described above, we find that strongly ordered states exist also in the original 
model. 

When T^ ax — > from below, the transient times increases in a much more erratic way 
than what was found for the transition T^ ax — > 1, and no simple behaviour was observed. 
(See Fig. 5 and 6). 

Much of the behavior described above can also be found in small systems, on which we 
can do some analytical calculations. Consider a system of only two sites, a and b. If T b is 
the state of b right after a has toppled, we can define the return map R b (T b ) as the state of 
b after the next toppling of a. For uniform thresholds it is easy to show that R b (T b ) = T b , 
where the "avalanche size" 5 = 1, i.e. there exist infinitely many marginally stable period 
two fixed points. However, if we introduce non-uniform thresholds by letting T* = 1 + e the 
return map will be 

R b (T b ) = T b + e(a + 1), (7) 

that is, the stress of b will increase with the amount e(a + 1) for each subsequent toppling of 
a and b. Of course, this can only continue until T b becomes larger than 1 — e(a + 1). Then, 
as b topples, T b — > and T a — >■ <p where a < <p < e(a + 1) + a. For a to topple next, the 
condition l + e- ^<lore<^ must be satisfied. A sufficient condition for this is that 
e < a, the necessary condition will depend on initial conditions. If the condition is satisfied, 
the following toppling of a will trigger b and we get an "avalanche" of size s = 2. We define 
the return map R a (<p) as the state of a after this avalanche and obtain 
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R a (<p) = a(l + e)(l + a)-a<p. (8) 

Since <p > a, R a (<p) < s(a + 1) + a and the next avalanche will also be of size s = 2. Thus, 
we have the general return map 

R n a +1 = a(l + £ )(1 + a) - aJC (9) 

This is a simple linear one- dimensional map with one fixed point, 

R a = a(l + e), (10) 

which is stable since a < 1. It is the same type of period one fixed point we found for the 
square lattices, but now there are only two A;- values, k = 0,1. From Eq.(7) we see that 

N(e,a)<x 1 (11) 
e(a + 1) 

where N is the number of avalanches before the fixed point is reached. The global increase 
of stress is proportional to the time, and from Eq.(7), 

A£ n _! + At n oc 1 + e - a, (12) 

where At n is the time between the n-th and the (n + l)-th avalanche. This yields a transient 
time 

<tr =Y At n GC — -. (13) 

When e — > this gives t tT ~ (r^ ax — l) -1 , which is in good agreement with the results 
of the square lattices. Still, we do not find it trivial that a lattice with local coupling has 
qualitatively the same asymptotic behavior as a system with only two sites. 

In summary, it is demonstrated that the coupled map lattice of Olami et al. is not robust 
as a model of earthquakes, since changing the threshold at just one site completely changes 
its behaviour. In a certain range of parameters, the resultant asymptotic states are highly 
organized, even with open boundary conditions. One may speculate if a similar doping of 
other degenerate type of models could have a similar effect. 

FT is grateful to Kim Christensen for several useful discussions. 
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FIG. 1. Distribution of stresses immediately after an avalanche in the period one i.e. asymptotic 
state for a 50 X 50 lattice with periodic boundary conditions, a = 0.20 and a = 0.04. 
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FIG. 2. Distribution of stresses after an avalanche in a quasi asymptotic state (after 1 X 10 7 
avalanches) for a 50 X 50 lattice with open boundary conditions, a = 0.20 and a = 0.04. 
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FIG. 3. Average size of the last 500 avalanches as a function of avalanche number for a 35 X 35 
lattice with periodic boundaries and one doped site, a = 0.1825 and T^ ax = 1.21. 
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FIG. 4. Average size of the last 500 avalanches as a function of avalanche number for a 35 X 35 
lattice with periodic boundaries and one doped site, a = 0.1825 and T^ ax = 1.40. 
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FIG. 5. Transient times for the strongly ordered state as a function of T^ ax — 1 for different 
values of a andi = 15. The dashed lines indicate the values T* aY = for the different a. 
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FIG. 6. Transient times for the strongly ordered state as a function of T^ ax — 1 for different 



values of L and a = 0.20. The dashed line indicates the value T* 



i 

l-a ' 



11 













P 




c 























10 20 30 40 50 

L 

FIG. 7. Transient times for a lattice with a = 0.22 and T^ ax = 1.01 as a function of lattice size 
L. For each L, the mean value for 20 different initial conditions has been plotted. The data points 
are well fitted to a straight line with slope ~ 1.6. 
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FIG. 8. Transient times as a function of a for L = 15 and T^ ax = 1.02. The initial conditions 
are the same for all values of a. 
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